Since release of the Biological Opinion on oil and gas activities in the Gulf of Mexico (NMFS 2020) that used a published density surface model (Roberts et al. 2016) to describe the distribution of the critically endangered Rice’s whale (Balaenoptera ricei), a new density surface model (Litz et al. 2022) has been made available. Importantly, this model extends the distribution of Rice’s whale beyond its initial core habitat in the Eastern Gulf of Mexico to the West, where it had previously only been acoustically detected (Soldevilla et al. 2022). This report replicates the Biological Opinion’s ship strike analysis using the newer Rice’s whale distributional model. Given the wider distribution of Rice’s whale, an alternative new Whale Area is suggested to reduce ship strike risk with the Rice’s whale based simply on location (25.5º N and higher) and depth (100 to 400 m).
2 Whale Densities
The new density surface model (Litz et al. 2022) uses approximately 40 km2 hexagons as its spatial unit to describe number of individuals per 40 km2 in a Lambert Conformal Conic projection, whereas the original model (Roberts et al. 2016) used 100 km2 cells in a custom equal area Albers projection to describe number of individuals per 100 km2. The spatial unit for this new analysis is also 100 km2 cells but in the web Mercator projection (EPSG:3857) in order to readily map results online with common “slippy” basemaps, like the Esri Ocean Basemap. All layers were clipped to the study area of the U.S. Exclusive Economic Zone (EEZ) within the Gulf of Mexico.
Normally converting polygons to raster extracts only the centroid point of the raster cell from the underlying polygon. In order to capture the entirety of the underlying geometric densities, a vector-based intersection was first performed on all layers (whale hexagons, ship cells, and new units) before summarizing to the raster cell as area-weighted means.
In order to adjust for slight differences from projecting coordinate reference systems and rounding errors, the new 100 km2 whale density grid was adjusted so the sum of individuals predicted throughout the study area is equal to 51.3, the most recent abundance estimate (Garrison, Ortega-Ortiz, and Rappucci 2020) from the same 2017 and 2018 surveys as the new density model (Litz et al. 2022) was derived.
Compared to the original distribution (Figure 1), the whales are now concentrated along the strip from 100 to 400 m extending into the Western Gulf of Mexico (Figure 2).
Figure 1: Map of previous whale densities (Roberts et al. 2016) as 100 km2 cells used by (NMFS 2020) showing the dominance in the northeastern corner of the Gulf of Mexico. The original Whale Area (p. 292 of NMFS 2020) is depicted by the pink outline polygon for vessel slowdown and nighttime avoidance. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
The original Whale Area is described (p. 292 of NMFS 2020) as:
This opinion defines the Bryde’s whale area to include the area from 100- to 400- meter isobaths from 87.5° W to 27.5° N as described in the status review (Rosel 2016) plus an additional 10 km around that area.
There was only a tiny marginal improvement in capturing additional whale densities by adding the 10 km buffer used to create the original Whale Area (pink outline in Figure 1). In generating the new Whale Area, the ease of navigation with simpler description in terms only of a southern limit and depth range outweighed this marginal improvement (red outline in Figure 2). This new Whale Area captures 94% of the population from the new density estimates (Litz et al. 2022) compared to only 52% of the original Whale Area (Table 1).
Figure 2: Map of new whale densities (Litz et al. 2022) as 100 km2 cells showing a distribution throughout the region. The newly recommended Whale Area is depicted by the red outline polygon for vessel slowdown and nighttime avoidance using similar logic as to (NMFS 2020). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Table 1:
Table of new whale densities (Litz et al. 2022) summarized by total
study area (U.S. Gulf of Mexico), previous Whale Area (NMFS 2020) and
newly proposed Whale Area.
Item
#
%
Whales in Study (U.S. Gulf of Mexico)
51
100%
Whales in Original Whale Area (NMFS, 2020)
27
52%
Whales in New Whale Area
48
94%
3 Vessel Traffic
In order to evaluate the threat of ship strike to Rice’s whales, we used the same AIS data from 2014 to 2018 as the Biological Opinion (NMFS 2020). This data is based on a grid of cells ~126 km2 in Albers equal area projection. Traffic in terms of kilometers (km) traversed within a cell was differentiated based on speed (≤ 10 knots or > 10 knots) and type (oil & gas or all types). In order to produce maps similar to the original Biological Opinion (NMFS 2020) showing spatial variation, colors were assigned to the Jenks natural breaks of the distribution of values (Figures 3, 4, 5, 6).
Figure 3: Map of annual average traffic (km) for all vessel types at all speeds from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Figure 4: Map of annual average traffic (km) for oil and gas vessels at all speeds from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Figure 5: Map of annual average traffic (km) for all vessel types > 10 knots from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Figure 6: Map of annual average traffic (km) for oil and gas vessels > 10 knots from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
4 Vessel Risk to Whales
The vessel risk (\(R\)) to whales is calculated as a simple multiplication of number of whales (\(W\)) and km of vessel traffic (\(V\)) (Equation 1).
\[
R = W * V
\tag{1}\]
This risk (\(R\)) can be further differentiated by vessel (\(V\)) type and speed (Figures 7, 8; Table Table 2).
Figure 7: Map of risk (# whales * km vessel traffic) for all vessels at all speeds. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Figure 8: Map of risk (# whales * km vessel traffic) for oil and gas vessels > 10 knots. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker).
Table 2:
Vessel strike risk (# whales * km vessel traffic) to Rice’s whales
for oil and gas vessels compared with all vessels.
Year
Vessel Strike Risk
All Vessels
Oil & Gas
%
All speeds
2015
92,849
38,428
41%
2016
85,323
33,280
39%
2017
86,236
33,858
39%
2018
100,326
38,668
39%
Avg
91,183
36,059
40%
> 10 knots
2015
71,621
23,153
32%
2016
67,878
21,519
32%
2017
68,494
21,909
32%
2018
79,759
25,209
32%
Avg
71,938
22,948
32%
Finally, we can evaluate the risk reduction of the original Whale Area proposed in the Biological Opinion (NMFS 2020) compared with the newly proposed Whale Area encompassing the wider Gulf of Mexico Table 3.
Table 3:
Reduction of vessel strike risk (# whales * km vessel traffic) to
Rice’s whales with enforcement of original (NMFS 2020) and new Whale
Areas. All percentage (%) reductions are compared to All Vessels for
given speeds.
Allaire, J. 2022. “Quarto: R Interface to’Quarto’Markdown Publishing System.”R Package Version 1.
Garrison, Lance, Joel Ortega-Ortiz, and Gina Rappucci. 2020. “Abundance of Marine Mammals in Waters of the U.S. Gulf of Mexico During the Summers of 2017 and 2018.” PRBD-2020-07.
Litz, Jenny, Laura Aichinger Dias, Gina Rappucci, Anthony Martinez, Melissa Soldevilla, Lance Garrison, and Keith Mullin. 2022. “Cetacean and Sea Turtle Spatial Density Model Outputs from Visual Observations Using Line-Transect Survey Methods Aboard NOAA Vessel and Aircraft Platforms in the Gulf of Mexico.”
Lowndes, Julia S. Stewart, Benjamin D. Best, Courtney Scarborough, Jamie C. Afflerbach, Melanie R. Frazier, Casey C. O’Hara, Ning Jiang, and Benjamin S. Halpern. 2017. “Our Path to Better Science in Less Time Using Open Data Science Tools.”Nature Ecology & Evolution 1 (6): 1–7. https://doi.org/10.1038/s41559-017-0160.
NMFS. 2020. “‘Biological Opinion on the Federally Regulated Oil and Gas Program Activities in the Gulf of Mexico,’ 13 March 2020, a Consultation Conducted by the Endangered Species Act Interagency Cooperation Division, Office of Protected Resources, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, U.S. Department of Commerce.”
Pebesma, Edzer. 2018. “The R Journal: Simple Features for R: Standardized Support for Spatial Vector Data.”The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
R Core Team, R. 2023. “R: A Language and Environment for Statistical Computing.” Vienna, Austria.
Roberts, Jason J., Benjamin D. Best, Laura Mannocci, Ei Fujioka, Patrick N. Halpin, Debra L. Palka, Lance P. Garrison, et al. 2016. “Habitat-Based Cetacean Density Models for the U.S. Atlantic and Gulf of Mexico.”Scientific Reports 6 (March): 22615. https://doi.org/10.1038/srep22615.
Soldevilla, Melissa S., Amanda J. Debich, Lance P. Garrison, John A. Hildebrand, and Sean M. Wiggins. 2022. “Rice’s Whales in the Northwestern Gulf of Mexico: Call Variation and Occurrence Beyond the Known Core Habitat.”Endangered Species Research 48 (July): 155–74. https://doi.org/10.3354/esr01196.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse.”Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Source Code
---title: "Spatial analysis of ship strike risk for Rice’s whale in the Gulf of Mexico"author: "Benjamin D. Best, Ph.D. (<ben@ecoquants.com>)"date: nowdate-format: "YYYY-MM-DD HH:mm (z)"bibliography: "ricei.bib"format: html: toc: true number-sections: true number-depth: 3 code-fold: true code-tools: true docx: toc: true toc-depth: 3 toc-title: "Contents" number-sections: true code-annotations: false execute: echo: false warning: falseeditor_options: chunk_output_type: console---```{r setup, include=FALSE}knitr::opts_chunk$set(echo = F,warning = F,message = F)source(here::here("scripts/functions.R"))source(here("scripts/paths.R"))ply_study <-read_sf(study_geo)ply_cells <-read_sf(cells_geo)ply_units <-read_sf(units_gpkg)ply_ships <-read_sf(ships_geo)ply_whales <-read_sf(whales_geo)ply_bia <-read_sf(bia_geo)ply_wab <-read_sf(wab_geo)ply_wan <-read_sf(wan_geo)tbl_ships <-read_csv(ships_csv)tbl_whales <-read_csv(whales_csv)lns_depth_contours <-read_sf(depth_contours_geo)stk_rast <-rast(rast_tif)r_cells <- stk_rast["cell_id"]tbl_units <- ply_units |>st_drop_geometry()``````{r}#| label: docx-fixes#| eval: !expr "!knitr::is_html_output()"#| results: asis# Manual docx: mv tbl outside figure table. Reference > Insert Caption. Add Table 1 and Insert Hyperlink to This Document "tbl-whales-by-area" to replace "?@tbl-whales-by-area".cat("## TODO: docx manual fixes {.unnumbered}- [ ] Replace all missed `gt` table references `?@` with `Table #` hyperlink.- [ ] Place cursor inside table, Reference > Insert Caption. Copy & paste from HTML output.- [ ] Drop `?(caption)` instances.- [ ] Move `gt` tables outside figure.")```## AbstractSince release of the Biological Opinion on oil and gas activities in the Gulf of Mexico [@nmfsBiologicalOpinionFederally2020] that used a published density surface model [@robertsHabitatbasedCetaceanDensity2016] to describe the distribution of the critically endangered Rice's whale (_Balaenoptera ricei_), a new density surface model [@litzCetaceanSeaTurtle2022] has been made available. Importantly, this model extends the distribution of Rice's whale beyond its initial core habitat in the Eastern Gulf of Mexico to the West, where it had previously only been acoustically detected [@soldevillaRiceWhalesNorthwestern2022]. This report replicates the Biological Opinion's ship strike analysis using the newer Rice's whale distributional model. Given the wider distribution of Rice's whale, an alternative new Whale Area is suggested to reduce ship strike risk with the Rice's whale based simply on location (25.5º N and higher) and depth (100 to 400 m).## Whale DensitiesThe new density surface model [@litzCetaceanSeaTurtle2022] uses approximately 40 km^2^ hexagons as its spatial unit to describe number of individuals per 40 km^2^ in a Lambert Conformal Conic projection, whereas the original model [@robertsHabitatbasedCetaceanDensity2016] used 100 km^2^ cells in a custom equal area Albers projection to describe number of individuals per 100 km^2^. The spatial unit for this new analysis is also 100 km^2^ cells but in the web Mercator projection (EPSG:3857) in order to readily map results online with common "slippy" basemaps, like the Esri Ocean Basemap. All layers were clipped to the study area of the U.S. Exclusive Economic Zone (EEZ) within the Gulf of Mexico.Normally converting polygons to raster extracts only the centroid point of the raster cell from the underlying polygon. In order to capture the entirety of the underlying geometric densities, a vector-based intersection was first performed on all layers (whale hexagons, ship cells, and new units) before summarizing to the raster cell as area-weighted means.In order to adjust for slight differences from projecting coordinate reference systems and rounding errors, the new 100 km^2^ whale density grid was adjusted so the sum of individuals predicted throughout the study area is equal to `r n_whales_Garrison2020`, the most recent abundance estimate [@garrisonAbundanceMarineMammals2020] from the same 2017 and 2018 surveys as the new density model [@litzCetaceanSeaTurtle2022] was derived. <!-- - discuss changing distribution -->Compared to the original distribution (@fig-map-whales-old), the whales are now concentrated along the strip from 100 to 400 m extending into the Western Gulf of Mexico (@fig-map-whales-new).```{r}#| label: fig-map-whales-old#| fig-cap: "Map of previous whale densities [@robertsHabitatbasedCetaceanDensity2016] as 100 km^2^ cells used by [@nmfsBiologicalOpinionFederally2020] showing the dominance in the northeastern corner of the Gulf of Mexico. The original Whale Area [p. 292 of @nmfsBiologicalOpinionFederally2020] is depicted by the pink outline polygon for vessel slowdown and nighttime avoidance. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast(r =rast(whales_roberts2016_img), legend_title ="Whales<br><small>(Roberts, 2016)<br>(# / 100 km<sup>2</sup>)</small>",group ="whales_per_100km2_Roberts2016",add_depth_contours =T,add_ply_wab = T)```The original Whale Area is described [p. 292 of @nmfsBiologicalOpinionFederally2020] as:> This opinion defines the Bryde’s whale area to include the area from 100- to 400- meter isobaths from 87.5° W to 27.5° N as described in the status review (Rosel 2016) plus an additional 10 km around that area.There was only a tiny marginal improvement in capturing additional whale densities by adding the 10 km buffer used to create the original Whale Area (pink outline in @fig-map-whales-old). In generating the new Whale Area, the ease of navigation with simpler description in terms only of a southern limit and depth range outweighed this marginal improvement (red outline in @fig-map-whales-new). This new Whale Area captures 94% of the population from the new density estimates [@litzCetaceanSeaTurtle2022] compared to only 52% of the original Whale Area (@tbl-whales-by-area).```{r}#| label: fig-map-whales-new#| fig-cap: "Map of new whale densities [@litzCetaceanSeaTurtle2022] as 100 km^2^ cells showing a distribution throughout the region. The newly recommended Whale Area is depicted by the red outline polygon for vessel slowdown and nighttime avoidance using similar logic as to [@nmfsBiologicalOpinionFederally2020]. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast(r = stk_rast["whales_avg_n_per_100km2"], legend_title ="Whales<br><small>(Litz, 2022)<br>(# / 100 km<sup>2</sup>)</small>",group ="whales_per_100km2_Litz2022",add_depth_contours=T,add_ply_wan = T)``````{r}#| label: tbl-whales-by-area#| tbl-cap: Table of new whale densities (Litz et al. 2022) summarized by total study area (U.S. Gulf of Mexico), previous Whale Area (NMFS 2020) and newly proposed Whale Area.a <-read_csv(whales_n_by_area_csv)d <-tribble(~Item, ~`#`, ~`%`,"Whales in Study (U.S. Gulf of Mexico)", a$n_whales_gom, 1,"Whales in Original Whale Area (NMFS, 2020)", a$n_whales_wab, a$pct_whales_wab,"Whales in New Whale Area", a$n_whales_wan, a$pct_whales_wan)d |>gt() |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent(`%`, decimals =0) |>opt_stylize(style =6)```## Vessel TrafficIn order to evaluate the threat of ship strike to Rice's whales, we used the same AIS data from 2014 to 2018 as the Biological Opinion [@nmfsBiologicalOpinionFederally2020]. This data is based on a grid of cells ~126 km^2^ in Albers equal area projection. Traffic in terms of kilometers (km) traversed within a cell was differentiated based on speed (≤ 10 knots or > 10 knots) and type (oil & gas or all types). In order to produce maps similar to the original Biological Opinion [@nmfsBiologicalOpinionFederally2020] showing spatial variation, colors were assigned to the Jenks natural breaks of the distribution of values (Figures [-@fig-ships-avg-all-gt01], [-@fig-ships-avg-boem-gt01], [-@fig-ships-avg-all-gt10], [-@fig-ships-avg-boem-gt10]).```{r}#| label: fig-ships-avg-all-gt01#| fig-cap: "Map of annual average traffic (km) for all vessel types at all speeds from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast_jenks(stk_rast["ships_avg_all_gt01"], "Traffic,<br><small>all types,<br>all speeds,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-boem-gt01#| fig-cap: "Map of annual average traffic (km) for oil and gas vessels at all speeds from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast_jenks(stk_rast["ships_avg_boem_gt01"], "Traffic,<br><small>oil & gas,<br>all speeds,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-all-gt10#| fig-cap: "Map of annual average traffic (km) for all vessel types > 10 knots from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast_jenks(stk_rast["ships_avg_all_gt01"], "Traffic,<br><small>all types,<br>> 10 knots,<br>avg year<br>(km)</small>", "RdYlGn")``````{r}#| label: fig-ships-avg-boem-gt10#| fig-cap: "Map of annual average traffic (km) for oil and gas vessels > 10 knots from AIS data (2014 to 2018). Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."map_rast_jenks(stk_rast["ships_avg_boem_gt01"], "Traffic,<br><small>oil & gas,<br>> 10 knots,<br>avg year<br>(km)</small>", "RdYlGn")```## Vessel Risk to WhalesThe vessel risk ($R$) to whales is calculated as a simple multiplication of number of whales ($W$) and km of vessel traffic ($V$) (@eq-risk).$$R = W * V$$ {#eq-risk}This risk ($R$) can be further differentiated by vessel ($V$) type and speed (Figures [-@fig-risk-avg-all-gt01], [-@fig-risk-avg-boem-gt10]; Table [@tbl-risk-overview]).```{r}#| label: fig-risk-avg-all-gt01#| fig-cap: "Map of risk (# whales * km vessel traffic) for all vessels at all speeds. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."stk_rast["risk_avg_all_gt01"] |>map_rast("Risk,<br><small>all types,<br>all speeds<br>(# * km)</small>", "RdYlGn")``````{r}#| label: fig-risk-avg-boem-gt10#| fig-cap: "Map of risk (# whales * km vessel traffic) for oil and gas vessels > 10 knots. Depth contours are shown in dash blacked lines for 100 m (finer) and 400 m (thicker)."stk_rast["risk_avg_boem_gt10"] |>map_rast("Risk,<br><small>oil & gas,<br>> 10 knots <br>(# * km)</small>", "RdYlGn")``````{r}#| label: tbl-risk-overview#| tbl-cap: "Vessel strike risk (# whales * km vessel traffic) to Rice's whales for oil and gas vessels compared with all vessels."d <-read_csv(risk_overview_csv) |>pivot_longer(-yr,names_to ="var",values_to ="val") |># table(d$var)separate_wider_delim( var,delim ="_",names =c("v","type","speed")) |>pivot_wider(names_from =c(v, type),names_sep ="_",values_from = val) |>relocate(speed) |>arrange(speed, yr)d <- d |>mutate(yr =as.character(yr)) |>bind_rows(bind_cols(yr ="Avg", d |>group_by(speed) |>summarise(across(where(is.numeric) &!yr, mean)) ) ) |>arrange(speed, yr)d <- d |>mutate(speed =case_match( speed, "gt01"~"All speeds","gt10"~"> 10 knots") )d |>group_by(speed) |>gt() |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent( pct_boem, decimals =0) |>opt_stylize(style =6) |>cols_label(yr ="Year",risk_all ="All Vessels",risk_boem ="Oil & Gas",pct_boem ="%") |>tab_spanner(label ="Vessel Strike Risk",columns =c( risk_all, risk_boem, pct_boem))```Finally, we can evaluate the risk reduction of the original Whale Area proposed in the Biological Opinion [@nmfsBiologicalOpinionFederally2020] compared with the newly proposed Whale Area encompassing the wider Gulf of Mexico [@tbl-risk-reduction-by-areas].```{r}#| label: tbl-risk-reduction-by-areas#| tbl-cap: "Reduction of vessel strike risk (# whales * km vessel traffic) to Rice's whales with enforcement of original (NMFS 2020) and new Whale Areas. All percentage (%) reductions are compared to All Vessels for given speeds."d <-read_csv(risk_reduction_by_areas_csv) |>select(yr, matches("wan|wab$")) |>rename_with(~str_replace(.x, "^pct_risk_reduced", "pctriskreduced"),starts_with("pct_risk_reduced")) |>pivot_longer(-yr,names_to ="var",values_to ="val") |># table(d$var)separate_wider_delim( var,delim ="_",names =c("v","type","speed","area")) |>pivot_wider(names_from =c(v, area),names_sep ="_",values_from = val) |>relocate(speed, type) |>arrange(speed, type, yr)# names(d) |> sort() |> paste(collapse="\n") |> cat()d <- d |>mutate(yr =as.character(yr)) |>bind_rows(bind_cols(yr ="Avg", d |>select(-yr) |>group_by(speed, type) |>summarise(across(where(is.numeric), mean),.groups ="drop") ) ) |>arrange(speed, type, yr)d <- d |>mutate(speed =case_match( speed, "gt01"~"All speeds","gt10"~"> 10 knots"),type =case_match( type, "all"~"All vessels","boem"~"Oil & Gas vessels") )d |>group_by(speed, type) |>gt(rowname_col ="year") |>fmt_number(decimals=0) |># c(risk_all, risk_boem), fmt_percent(c( pctriskreduced_wab, pctriskreduced_wan), decimals =0) |>opt_stylize(style =6) |>cols_label(yr ="Year",risk_wab ="Original",pctriskreduced_wab ="%",risk_wan ="New",pctriskreduced_wan ="%") |>tab_spanner(label ="Risk Reduction by Area",columns =c( risk_wab, pctriskreduced_wab, risk_wan, pctriskreduced_wan)) |>cols_align(align ="right",columns = yr)```## Reproducible ResultsThis report was produced using the principles of reproducible research [@lowndesOurPathBetter2017] with the R programming language [@rcoreteamLanguageEnvironmentStatistical2023]. Statistical analysis were performed using the libraries and methods of the [`tidyverse`](https://www.tidyverse.org/)[@wickhamWelcomeTidyverse2019] and spatial features [`sf`](https://r-spatial.github.io/sf/index.html)[@pebesmaJournalSimpleFeatures2018] output to a [Quarto](https://quarto.org/) document [@allaireQuartoInterfaceQuarto2022]. All source code is available in the Github repository [github.com/ecoquants/ricei](https://github.com/ecoquants/ricei). The interactive version of this report is available at [ecoquants.com/ricei](https://ecoquants.com/ricei).## References {.unnumbered}